---
title: "Effects of environmental stressors on daily governance"
subtitle: Main text replication code
author: Nick Obradovich, Dustin Tingley, and Iyad Rahwan
output: pdf_document
---

<!--
Rscript -e "rmarkdown::render('./main_text_replication_code.Rmd')"
-->

```{r opts, echo=FALSE, message=FALSE, warning=FALSE, include=FALSE}
knitr::opts_chunk$set(fig.path='.', warning=F, message=F, fig.retina=NULL, cache=T, cache.lazy=F, autodep=T, echo=F, eval=T)
```

```{r packages}
# packages----
library(bit64) # long integers
library(ggthemes) #themes for plotting
library(sandwich) #robust se estimator
library(lmtest) #lm functions
library(ordinal) #ologit funcitons
library(memisc) #import stata data
library(ggmap) #make nice plots with ggplot mapping
library(ggplot2) #nice plots
library(Hmisc) #multipurpose package
library(RColorBrewer) #make nice colors
library(grid) #plotting package
library(MASS) #multipurpose package
library(gplots) #nice plotting
library(lfe) #run linear fixed effects models
library(stargazer) #output nice tables
library(lubridate) #deal with dates
library(zoo) #deal with dates
library(reshape2) #powerful reshape package
library(sp) #general spatial class package
library(gstat) #geostat package
library(spacetime) #store data in proper spatial/temporal class
library(raster) #convert to raster and vice versa
library(maptools) #plot and modify shapefiles
library(rgeos) #for use with maptools for shapefiles
library(rgdal) #for reading in shapefiles
library(RSAGA) #needed for spatial downscaling
library(RCurl) #also needed for spatial downscaling
library(dplyr) #for data frame manipulation functions
library(plyr) #for additional manipulation functions
library(data.table) #for data.table functions
library(foreign) #for reading in .dta files
library(parallel) #for parallel processing functions
library(foreach) #for plyr parallel
library(doMC) #run multicore processes
library(caTools) #functions for fast running mean and sd
library(ncdf4) #tools for netCDF packages
library(stringr) #string tools
library(knitr) #knitting
library(pander) #pandering
library(rasterVis) # plotting rasters
library(GGally) # for correlation plots
library(geosphere) # spherical distance
library(Rcpp) # for C++
library(RcppArmadillo) # for C++ (and Conley errors)
library(matrixStats) # matrix calculations
library(ggExtra) # for marginal histograms
library(gridExtra) # for multiplots
library(viridis) # for colors
library(Cairo) # pdf device
library(gtable) # for gtable filter
library(devtools)
#devtools::install_github("wilkelab/cowplot")
library(cowplot)
library(captioner)

# set cores for parallel processing----
registerDoMC(10)
```

```{r functions}
# functions for grabbing p-values and coefficients from summaries of felm models----
coef <- function(reg, name) {
    if (round(reg$coefficients[rownames(reg$coefficients)==name, 1], 3) == 0) {
        format(round(reg$coefficients[rownames(reg$coefficients)==name, 1], 4), scientific=FALSE) 
    }
    else round(reg$coefficients[rownames(reg$coefficients)==name, 1], 3)
}
pval <- function(reg, name) {
  if (round(reg$coefficients[rownames(reg$coefficients)==name, 4], 3) >=0.001) 
  {round(reg$coefficients[rownames(reg$coefficients)==name, 4], 3)} 
  else paste("< 0.001")
}

# captioner setting----
fig_num <- captioner(prefix="Fig.")
tab_num <- captioner(prefix="Table")
eq_num <- captioner(prefix="Equation")
sfig_num <- captioner(prefix="Fig. S", auto_space=FALSE)
stab_num <- captioner(prefix="Table S", auto_space=FALSE)
seq_num <- captioner(prefix="Equation S", auto_space=FALSE)
```

```{r data}
# public data----
load("./replication.RData")

# proprietary data----
# d <- readRDS('./hazel_weather_inspection_balanced_panel.rds.gz')
# hazel <- readRDS('./hazel_geolocated_master.rds.gz')
# hazel_county <- readRDS('./hazel_counties.rds.gz')
```

```{r models}
# police stops----
pol <- summary(readRDS("./log_num_stop_felm.rds.gz"))

# fatal crashes----
cr <- summary(readRDS("./prob_crash_felm.rds.gz"))

# food safety inspections----
insp <- readRDS("./facility_day_felm_quadratic_summary.rds.gz")

# food safety violations----
viol <- summary(readRDS("./num_violation_felm_quad.rds.gz"))
```

# Figure 1

```{r fig1, echo=TRUE, eval=FALSE}
# plot time-series of main data----
## first calculate relevant time series
stoptime <- stops[, .(tot_stops = sum(num_stop, na.rm=T)), by=.(wkyr = as.Date(cut(date, breaks = "week")))]
crashtime <- crash[date < as.Date("2015-12-28"), .(num_crash = sum(accidents, na.rm=T)), by=.(wkyr = as.Date(cut(date, breaks = "week")))]
insptime <- hazel[date < as.Date("2016-12-01"), .(tot_inspections = sum(uniqueN(inspection_id), na.rm=T)), by=.(wkyr = as.Date(cut(date, breaks = "week")))]
# violtime <- hazel[, .(violperinsp = sum(!is.na(violation_id), na.rm=T)), by=.(inspection_id, date)]
# violtime <- violtime[date < as.Date("2016-12-01"), .(violperinsp = sum(violperinsp, na.rm=T)), by=.(wkyr = as.Date(cut(date, breaks = "week")))]
violtime <- hazel[date < as.Date("2016-12-01"), .(tot_violations = sum(uniqueN(violation_id), na.rm=T)), by=.(wkyr = as.Date(cut(date, breaks = "week")))]

## then plot
stop.plot <- ggplot(stoptime, aes(x=wkyr, y=tot_stops)) +
                    ylab("Police Stops") +
                    xlab("Date") +
                    geom_line(color='gray50', size=0.35) + 
                    geom_area(fill="blue4", alpha=0.25) +
                    theme(plot.title = element_text(face="bold", size=40),
                          axis.title.x = element_text(size=20, face="bold"),
                          axis.title.y = element_text(size=20, face="bold"),
                          axis.text.y = element_text(size=15, face="bold"),
                          axis.text.x = element_text(size=15, face="bold"),
                          axis.ticks.y = element_line(),
                          panel.grid.major = element_blank(),
                          panel.grid.minor = element_blank(),
                          panel.background = element_blank(),
                          axis.line.x = element_line(),
                          axis.line.y = element_line()
                    )

crash.plot <- ggplot(crashtime, aes(x=wkyr, y=num_crash)) +
                    ylab("Fatal Crashes") +
                    xlab("Date") +
                    geom_line(color='gray50', size=0.35) + 
                    geom_area(fill="red4", alpha=0.25) +
                    scale_y_continuous(labels = function(x) formatC(x, format = "e", digits = 0)) +
                    theme(plot.title = element_text(face="bold", size=40),
                          axis.title.x = element_text(size=20, face="bold"),
                          axis.title.y = element_text(size=20, face="bold"),
                          axis.text.y = element_text(size=15, face="bold"),
                          axis.text.x = element_text(size=15, face="bold"),
                          axis.ticks.y = element_line(),
                          panel.grid.major = element_blank(),
                          panel.grid.minor = element_blank(),
                          panel.background = element_blank(),
                          axis.line.x = element_line(),
                          axis.line.y = element_line()
                    )

insp.plot <- ggplot(insptime, aes(x=wkyr, y=tot_inspections)) +
                    ylab("Food Safety Inspections") +
                    xlab("Date") +
                    geom_line(color='gray50', size=0.35) + 
                    geom_area(fill="orange", alpha=0.25) +
                    scale_y_continuous(labels = function(x) format(x, scientific = TRUE)) +
                    theme(plot.title = element_text(face="bold", size=40),
                          axis.title.x = element_text(size=20, face="bold"),
                          axis.title.y = element_text(size=20, face="bold"),
                          axis.text.y = element_text(size=15, face="bold"),
                          axis.text.x = element_text(size=15, face="bold"),
                          axis.ticks.y = element_line(),
                          panel.grid.major = element_blank(),
                          panel.grid.minor = element_blank(),
                          panel.background = element_blank(),
                          axis.line.x = element_line(),
                          axis.line.y = element_line()
                    )

viol.plot <- ggplot(violtime, aes(x=wkyr, y=tot_violations)) +
                    ylab("Food Safety Violations") +
                    xlab("Date") +
                    geom_line(color='gray50', size=0.35) + 
                    geom_area(fill="green4", alpha=0.25) +
                    scale_y_continuous(labels = function(x) formatC(x, format = "e", digits = 0)) +
                    theme(plot.title = element_text(face="bold", size=40),
                          axis.title.x = element_text(size=20, face="bold"),
                          axis.title.y = element_text(size=20, face="bold"),
                          axis.text.y = element_text(size=15, face="bold"),
                          axis.text.x = element_text(size=15, face="bold"),
                          axis.ticks.y = element_line(),
                          panel.grid.major = element_blank(),
                          panel.grid.minor = element_blank(),
                          panel.background = element_blank(),
                          axis.line.x = element_line(),
                          axis.line.y = element_line()
                    )

# maps of data----
## police stops map
names(counties)[6] <- "county_fips" # fix names
counties$county_fips <- as.character(counties$county_fips)
stops[, county_fips := as.character(county_fips)]
stop.county <- stops[, .(log_num_stop = log(sum(num_stop, na.rm=T))), by=.(county_fips, state)]
stop.county[is.infinite(log_num_stop), log_num_stop := 0]
stop.county <- sp::merge(counties, stop.county, by="county_fips")

p1 <- spplot(
       obj=stop.county,
       zcol=c("log_num_stop"),
       strip=FALSE,
       col.regions=colorRampPalette(c("gray80","blue4"))(100), # set colors
       as.table = TRUE,
       col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2)), # make legend on bottom of plot
       at=seq(min(stop.county$log_num_stop, na.rm=TRUE), max(stop.county$log_num_stop, na.rm=TRUE)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       sub=list("ln(Police Stops)", cex=1.4)
       )
p1 <- update(p1, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p1 <- p1 + layer_(sp.polygons(obj=us,
                              fill=gray(0.95),
                              col=gray(0.85)))

## crash map
names(counties)[6] <- "county_fips" # fix names
crash.county <- crash[, .(log_num_crash = log(sum(accidents, na.rm=T))), by=.(county_fips)]
crash.county[is.infinite(log_num_crash), log_num_crash := 0]
crash.county <- sp::merge(counties, crash.county, by="county_fips")

p2 <- spplot(
       obj=crash.county,
       zcol=c("log_num_crash"),
       strip=FALSE,
       col.regions=colorRampPalette(c("gray80","red4"))(100),
       as.table = TRUE,
       col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2)), # make legend on bottom of plot
       at=seq(min(crash.county$log_num_crash, na.rm=TRUE), max(crash.county$log_num_crash, na.rm=TRUE)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       sub=list("ln(Fatal Crashes)", cex=1.4)
       )
p2 <- update(p2, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p2 <- p2 + layer_(sp.polygons(obj=us,
                              fill=gray(0.95),
                              col=gray(0.85)))

## inspections map
names(counties)[6] <- "county_fips" # fix names
points <- hazel[, .(num_inspections = uniqueN(inspection_id)), by=.(facility_id, lat, lon)]
coordinates(points) <- ~lon+lat
projection(points) <- projection(counties)
points <- cbind(as.data.table(points), as.data.table(over(points, counties))$county_fips)
setnames(points, "V2","county_fips")
insp.county <- points[, .(log_num_insp = log(sum(num_inspections, na.rm=T))), by=.(county_fips)]
insp.county <- sp::merge(counties, insp.county, by="county_fips")

p3 <- spplot(
       obj=insp.county,
       zcol=c("log_num_insp"),
       strip=FALSE,
       col.regions=colorRampPalette(c("gray80","orange","orange4"))(100),
       as.table = TRUE,
       col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2)), # make legend on bottom of plot
       at=seq(min(insp.county$log_num_insp, na.rm=TRUE), max(insp.county$log_num_insp, na.rm=TRUE)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       sub=list("ln(Food Safety Inspections)", cex=1.4)
       )
p3 <- update(p3, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p3 <- p3 + layer_(sp.polygons(obj=us,
                              fill=gray(0.95),
                              col=gray(0.85)))

## violations map
names(counties)[6] <- "county_fips" # fix names
points <- hazel[, .(num_violations = uniqueN(violation_id)), by=.(facility_id, lat, lon)]
coordinates(points) <- ~lon+lat
projection(points) <- projection(counties)
points <- cbind(as.data.table(points), as.data.table(over(points, counties))$county_fips)
setnames(points, "V2","county_fips")
viol.county <- points[, .(num_violations = sum(num_violations, na.rm=T)), by=.(county_fips)]
viol.county[num_violations==0, num_violations := NA] # remove counties that don't record violations
viol.county[, log_num_viol := log(num_violations)]
viol.county <- sp::merge(counties, viol.county, by="county_fips")

p4 <- spplot(
       obj=viol.county,
       zcol=c("log_num_viol"),
       strip=FALSE,
       col.regions=colorRampPalette(c("gray80","green4"))(100), # set colors
       as.table = TRUE,
       col=NA,
       xlab=NULL, # kill x-axis
       ylab=NULL, # kill y-axis
       scales=list(draw=FALSE),
       margin=FALSE, # remove margin padding
       colorkey=list(space="bottom", labels=list(cex=1.2)), # make legend on bottom of plot
       at=seq(min(insp.county$log_num_insp, na.rm=TRUE), max(insp.county$log_num_insp, na.rm=TRUE)+0.000001, length.out=100),
       par.settings = plasmaTheme(),
       sub=list("ln(Food Safety Violations)", cex=1.4)
       )
p4 <- update(p4, par.settings = list(axis.line = list(col = "transparent"),
                                   par.main.text=list(cex=2),
                                   panel.background=list(col="transparent")))
p4 <- p4 + layer_(sp.polygons(obj=us,
                              fill=gray(0.95),
                              col=gray(0.85)))
# plot----
cairo_pdf(filename = "./figure1.pdf", width=18.33, height=5.83, fallback_resolution=1200)
grid.arrange(grobs = list(
                            arrangeGrob(rectGrob(gp=gpar(col="white")), p1, 
                                        rectGrob(gp=gpar(col="white")), p2, 
                                        rectGrob(gp=gpar(col="white")), p3, 
                                        rectGrob(gp=gpar(col="white")), p4, 
                                        widths=c(0.04, 0.21, 0.04, 0.21, 0.04, 0.21, 0.04, 0.21)),
                            arrangeGrob(stop.plot, crash.plot, insp.plot, viol.plot, ncol=4, nrow=1)
                            ), 
             ncol=1, nrow=2, heights=c(0.45, 0.55))
grid.text("a", x=unit(0.04, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("b", x=unit(0.29, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("c", x=unit(0.54, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("d", x=unit(0.79, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
dev.off()
```

```{r maps, out.width="100%", fig.cap="Spatial and temporal variation. This figure depicts the US counties covered by each of our data sources as well as the weekly variation in each series. Panel (a) plots the county locations of the over 70 million police stops in our data and illustrates that our police stops occur mostly between 2010 and 2015. Panel (b) displays the full US coverage of the more than 500,000 fatal crashes in our data and illustrates the seasonal trends in fatal crashes. Panel (c) shows the county distribution of the over 750,000 unique facilities across over four million inspections in our food safety data between 2012 and 2016. Panel (d) shows the same county distribution as panel (c) and highlights the fact that there are nearly three food safety violations per inspection on average in our data."}
knitr::include_graphics("figure1.pdf")
```

```{r}
maps <- fig_num(name = "maps", caption = "")
```

# Figure 2

```{r fig2, echo=TRUE, eval=FALSE}
# tmax----
t <- stops[tmax >= -15 & tmax <=45, seq(from=min(tmax, na.rm=T), to=max(tmax, na.rm=T), by=0.001)]
insp.t <- pol$coef[1]*t + pol$coef[2]*t^2
cr.t <- cr$coef[1]*t + cr$coef[2]*t^2 + 0.12 - max(insp.t) # slight re-shift for tangency
insp.t <- insp.t - max(insp.t) # put max of parabola at zero
tmax <- stops[, sample(tmax, length(t))]
t.data <- as.data.table(cbind(t, insp.t, cr.t, tmax))

t.plot <- ggplot(t.data) + 
    ylab("ln(Num. Police Stops) \n and P(Fatal Crash)") +
    xlab("Max. Temperature in \u2103") +
    geom_ribbon(aes(x=t, ymin=insp.t, ymax=cr.t, fill="gray85")) +
    geom_hline(yintercept = 0, linetype = 2, size = 0.1) +
    geom_line(aes(x=t, y=insp.t, color="cornflowerblue"), size=3, alpha=0.75) + 
    geom_line(aes(x=t, y=cr.t, color="red4"), size=3, alpha=0.75) + 
    coord_cartesian(xlim=c(min(t), max(t)), ylim=c(-0.75, 1.3), expand=TRUE) +
    scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
    scale_color_manual(limits=c("cornflowerblue", "red4"), values=c("cornflowerblue", "red4"), labels=c("Police Stops","P(Fatal Crash)")) +
    scale_fill_manual(limits=c("gray85"), values=c("gray85"), labels=c("Marginal Regulatory Gap")) +
    guides(color = guide_legend(override.aes = list(fill = NULL), order=1), fill = guide_legend(order=2)) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.margin = margin(-0.6,0,0,0, unit="cm"),
          legend.position = c(.25, 0.75),
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )
hist <- axis_canvas(t.plot, axis = "x") + 
        geom_histogram(data = t.data, aes(x = tmax), bins=20, fill='red4', color='gray60', alpha=0.75, size=0.1)
t.plot.full <- insert_xaxis_grob(t.plot, hist, position = "bottom", height = grid::unit(0.1, "null"))
t.plot.full <- ggdraw(t.plot.full)

# justy----
test <- ggplot(t.data) + 
    ylab("ln(Num. Police Stops) \n and P(Fatal Crash) in Percentage Points") +
    xlab("Max. Temperature in \u2103") +
    geom_ribbon(aes(x=t, ymin=insp.t, ymax=cr.t, fill="gray85")) +
    geom_hline(yintercept = 0, linetype = 2, size = 0.1) +
    geom_line(aes(x=t, y=insp.t, color="cornflowerblue"), size=3, alpha=0.75) + 
    geom_line(aes(x=t, y=cr.t, color="red4"), size=3, alpha=0.75) + 
    coord_cartesian(xlim=c(min(t), max(t)), ylim=c(-0.75, 1.3), expand=TRUE) +
    scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
    scale_color_manual(limits=c("cornflowerblue", "red4"), values=c("cornflowerblue", "red4"), labels=c("Police Stops","P(Fatal Crash)")) +
    scale_fill_manual(limits=c("gray85"), values=c("gray85"), labels=c("Regulatory Gap")) +
    guides(color = guide_legend(override.aes = list(fill = NULL), order=1), fill = guide_legend(order=2)) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_text(size=22, angle=90, vjust=-0.5, face="bold"),
          axis.text.y = element_text(size=22, angle=0, face="bold"),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_line(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          legend.title=element_blank(),
          legend.position = c(.75, 0.075),
          legend.text=element_text(size=16, face="bold"),
          legend.key = element_rect(colour = "black")
    )
hist <- axis_canvas(test, axis = "x") + 
        geom_histogram(data = t.data, aes(x = tmax), bins=20, fill='red4', color='gray60', alpha=0.75, size=0.1)
test <- insert_xaxis_grob(test, hist, position = "bottom", height = grid::unit(0.1, "null"))
justy <- gtable_filter(test, 'axis-l|ylab', trim=F)
justy <- ggdraw(justy)

# prcp----
p <- stops[, seq(from=min(prcp, na.rm=T), to=15, by=0.001)]
insp.p <- pol$coef[3]*p
cr.p <- cr$coef[3]*p # slight re-shift for tangency
insp.p <- insp.p - max(insp.p) # put max of parabola at zero
prcp <- stops[, sample(prcp, length(p))]
p.data <- as.data.table(cbind(p, insp.p, cr.p, prcp))

p.plot <- ggplot(p.data) + 
    ylab("ln(Num. Police Stops) \n and P(Fatal Crash)") +
    xlab("Precipitation in mm.") +
    geom_ribbon(aes(x=p, ymin=insp.p, ymax=cr.p, fill="gray85")) +
    geom_hline(yintercept = 0, linetype = 2, size = 0.1) +
    geom_line(aes(x=p, y=insp.p, color="cornflowerblue"), size=3, alpha=0.75) + 
    geom_line(aes(x=p, y=cr.p, color="red4"), size=3, alpha=0.75) + 
    coord_cartesian(xlim=c(min(p), 15), ylim=c(-0.75, 1.3), expand=TRUE) +
    scale_x_continuous(breaks=c(0, 5, 10, 15)) +
    scale_color_manual(limits=c("cornflowerblue", "red4"), values=c("cornflowerblue", "red4"), labels=c("Police Stops","P(Fatal Crash)")) +
    scale_fill_manual(limits=c("gray85"), values=c("gray85"), labels=c("Regulatory Gap")) +
    guides(color = guide_legend(override.aes = list(fill = NULL))) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.margin = margin(-0.6,0,0,0, unit="cm"),
          legend.position = "none",
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )
hist <- axis_canvas(p.plot, axis = "x") + 
        geom_histogram(data = p.data, aes(x = prcp), bins=16, fill='red4', color='gray60', alpha=0.75, size=0.1)
p.plot.full <- insert_xaxis_grob(p.plot, hist, position = "bottom", height = grid::unit(0.1, "null"))
p.plot.full <- ggdraw(p.plot.full)

# plot----
cairo_pdf(filename = "./figure2.pdf", width=15, height=7.5, fallback_resolution=1200)
grid.arrange(grobs=list(
  justy,
  arrangeGrob(t.plot.full, p.plot.full, ncol=2)),
  widths=c(0.08, 0.92), ncol=2) # plot all together
grid.text("a", x=unit(0.10, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("b", x=unit(0.56, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
dev.off()
```

```{r stops, out.width="100%", fig.cap="Marginal effects of temperature and precipitation on number of police stops and fatal vehicular crashes. Panel (a) depicts the relationship between county-level daily maximum temperatures and the log daily number of police stops in blue and between county-level daily maximum temperatures and probability of fatal vehicular crashes in percentage points in red. We draw from the estimation of Equation 1 and plot the predicted change in each outcome across the maximum temperature range within the data. Panel (b) also draws on the estimation of Equation 1 and depicts the effect of daily precipitation levels on our outcome variables. Due to the fixed effects in our analysis, the placement of these marginal effects curves on the y-axis is arbitrary; we set the curves tangent to one another for ease of inspection of the divergence of their slopes. Gray shaded areas depict the gap in marginal regulatory effort created by the divergence of these relationships under more extreme meteorological conditions. See SI: Marginal effects and SI: Regression tables for depiction of the 95\\% confidence intervals associated with the marginal effects of these estimates."}
knitr::include_graphics("figure2.pdf")
```

```{r}
stop <- fig_num(name = "stop", caption = "")
```

# Figure 3

```{r fig3, echo=TRUE, eval=FALSE}
# tmax ----
t <- d[tmax >= -15 & tmax <=45, seq(from=min(tmax, na.rm=T), to=max(tmax, na.rm=T), by=0.001)]
insp.t <- insp$coef[1]*t + insp$coef[2]*t^2
viol.t <- viol$coef[1]*t + viol$coef[2]*t^2 + 0.01 - max(insp.t) # slight re-shift for tangency
insp.t <- insp.t - max(insp.t) # put max of parabola at zero
tmax <- d[, sample(tmax, length(t))]
t.data <- as.data.table(cbind(t, insp.t, viol.t, tmax))

t.plot <- ggplot(t.data) + 
    ylab("P(Inspection) in Percentage Points\n and Number of Violations per Inspection") +
    xlab("Max. Temperature in \u2103") +
    geom_ribbon(aes(x=t, ymin=insp.t, ymax=viol.t, fill="gray85")) +
    geom_hline(yintercept = 0, linetype = 2, size = 0.1) +
    geom_line(aes(x=t, y=insp.t, color="orange"), size=3, alpha=0.75) + 
    geom_line(aes(x=t, y=viol.t, color="green4"), size=3, alpha=0.75) + 
    coord_cartesian(xlim=c(min(t), max(t)), ylim=c(-0.09, 0.16), expand=TRUE) +
    scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
    scale_color_manual(limits=c("orange", "green4"), values=c("orange", "green4"), labels=c("P(Inspection)", "Violations")) +
    scale_fill_manual(limits=c("gray85"), values=c("gray85"), labels=c("Marginal Regulatory Gap")) +
    guides(color = guide_legend(override.aes = list(fill = NULL), order=1), fill = guide_legend(order=2)) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.margin = margin(-0.6,0,0,0, unit="cm"),
          legend.position = c(.25, 0.75),
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )
hist <- axis_canvas(t.plot, axis = "x") + 
        geom_histogram(data = t.data, aes(x = tmax), bins=20, fill='orange', color='gray60', alpha=0.75, size=0.1)
t.plot.full <- insert_xaxis_grob(t.plot, hist, position = "bottom", height = grid::unit(0.1, "null"))
t.plot.full <- ggdraw(t.plot.full)

# justy----
test <- ggplot(t.data) + 
    ylab("P(Inspection) in Percentage Points\n and Number of Violations per Inspection") +
    xlab("Max. Temperature in \u2103") +
    geom_hline(yintercept = 0, linetype = 2, size = 0.1) +
    geom_line(aes(x=t, y=insp.t, color="orange"), size=3, alpha=0.75) + 
    geom_line(aes(x=t, y=viol.t, color="green4"), size=3, alpha=0.75) + 
    geom_ribbon(aes(x=t, ymin=insp.t, ymax=viol.t), alpha=0.1) +
    coord_cartesian(xlim=c(min(t), max(t)), ylim=c(-0.09, 0.16), expand=TRUE) +
    scale_x_continuous(breaks=c(-10, 0, 10, 20, 30, 40)) +
    scale_color_manual(limits=c("orange", "green4"), values=c("orange", "green4"), labels=c("P(Inspection)", "Violations")) +
    guides(color = guide_legend(override.aes = list(fill = NULL))) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_text(size=22, angle=90, vjust=-0.5, face="bold"),
          axis.text.y = element_text(size=22, angle=0, face="bold"),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_line(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          legend.title=element_blank(),
          legend.position = c(.75, 0.075),
          legend.text=element_text(size=16, face="bold"),
          legend.key = element_rect(colour = "black")
    )
hist <- axis_canvas(test, axis = "x") + 
        geom_histogram(data = t.data, aes(x = tmax), bins=20, fill='orange', color='gray60', alpha=0.75, size=0.1)
test <- insert_xaxis_grob(test, hist, position = "bottom", height = grid::unit(0.1, "null"))
justy <- gtable_filter(test, 'axis-l|ylab', trim=F)
justy <- ggdraw(justy)

# prcp----
p <- d[, seq(from=min(prcp, na.rm=T), to=15, by=0.001)]
insp.p <- insp$coef[3]*p
viol.p <- viol$coef[2]*p # slight re-shift for tangency
prcp <- d[, sample(prcp, length(p))]
p.data <- as.data.table(cbind(p, insp.p, viol.p, prcp))

p.plot <- ggplot(p.data) + 
    ylab("P(Inspection) in Percentage Points\n and Number of Violations per Inspection") +
    xlab("Precipitation in mm.") +
    geom_ribbon(aes(x=p, ymin=insp.p, ymax=viol.p, fill="gray85")) +
    geom_hline(yintercept = 0, linetype = 2, size = 0.1) +
    geom_line(aes(x=p, y=insp.p, color="orange"), size=3, alpha=0.75) + 
    geom_line(aes(x=p, y=viol.p, color="green4"), size=3, alpha=0.75) + 
    coord_cartesian(xlim=c(min(p), 15), ylim=c(-0.09, 0.16), expand=TRUE) +
    scale_x_continuous(breaks=c(0, 5, 10, 15)) +
    scale_color_manual(limits=c("orange", "green4"), values=c("orange", "green4"), labels=c("P(Inspection)", "Violations")) +
    scale_fill_manual(limits=c("gray85"), values=c("gray85"), labels=c("Regulatory Gap")) +
    guides(color = guide_legend(override.aes = list(fill = NULL))) +
    theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=25, angle=0, face="bold"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.text.x = element_text(size=20, face="bold"),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_blank(),
          legend.title=element_blank(),
          legend.margin = margin(-0.6,0,0,0, unit="cm"),
          legend.position = "none",
          legend.text=element_text(size=15, face="bold"),
          legend.key = element_rect(colour = "black")
    )
hist <- axis_canvas(p.plot, axis = "x") + 
        geom_histogram(data = p.data, aes(x = prcp), bins=16, fill='orange', color='gray60', alpha=0.75, size=0.1)
p.plot.full <- insert_xaxis_grob(p.plot, hist, position = "bottom", height = grid::unit(0.1, "null"))
p.plot.full <- ggdraw(p.plot.full)

# plot----
cairo_pdf(filename = "./figure3.pdf", width=15, height=7.5, fallback_resolution=1200)
grid.arrange(grobs=list(
  justy,
  arrangeGrob(t.plot.full, p.plot.full, ncol=2)),
  widths=c(0.08, 0.92), ncol=2) # plot all together
grid.text("a", x=unit(0.10, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
grid.text("b", x=unit(0.56, "npc"), y=unit(0.97, "npc"), gp=gpar(fontsize=33,font=2))
dev.off()
```

```{r insp, out.width="100%", fig.cap="Marginal effects of temperature and precipitation on food safety inspections and violations. Panel (a) depicts the relationship between daily maximum temperatures and the percentage point probability of facility-level food safety inspections in gold and between daily maximum temperatures and the number of food safety violations per facility-level inspection in green. It draws from the estimation of Equation 1 and plots the predicted change in each outcome across the maximum temperature range within the data. Panel (b) also draws on the estimation of Equation 1 and depicts the effect of daily precipitation levels on our outcome variables. Due to the fixed effects in our analysis, the placement of these marginal effects curves on the y-axis is arbitrary; we set the curves tangent to one another for ease of inspection of the divergence of their slopes. Gray shaded areas depict the gap in marginal regulatory effort created by the divergence of these relationships. See SI: Marginal effects and SI: Regression tables for the 95\\% confidence intervals associated with the marginal effects of these estimates."}
knitr::include_graphics("figure3.pdf")
```

```{r}
insps <- fig_num(name = "insps", caption = "")
```

# Figure 4

```{r fig4, echo=TRUE, eval=FALSE}
# calculate RCP8.5 differences for plotting----
id <- c('2050','2099') # set year id

## stops
meandiff.pol <- stack(mean(100*(exp(1)^diff2050.pol - 1)), mean(100*(exp(1)^diff2099.pol - 1))) # change in mean daily police stops in percentage points (convert changes in fitted logs to percentage points)
meandiff.pol <- setZ(meandiff.pol, id) # set year as z variable
names(meandiff.pol) <- c('2050','2099') # set names
meandiff.pol <- trim(meandiff.pol, padding=0, values=NA) # trim outer NAs

## crashes
meandiff.cr <- stack(mean(diff2050.cr), mean(diff2099.cr)) # change in mean daily fatal crash risk in percentage points
meandiff.cr <- setZ(meandiff.cr, id) # set year as z variable
names(meandiff.cr) <- c('2050','2099') # set names
meandiff.cr <- trim(meandiff.cr, padding=0, values=NA) # trim outer NAs

## food inspections
meandiff.insp <- stack(mean(diff2050.insp*1000), mean(diff2099.insp*1000)) # change in mean daily inspections per 100k facilities (multiply percentage points by 1000)
meandiff.insp <- setZ(meandiff.insp, id) # set year as z variable
names(meandiff.insp) <- c('2050','2099') # set names
meandiff.insp <- trim(meandiff.insp, padding=0, values=NA) # trim outer NAs

##food violations
meandiff.viol <- stack(mean(diff2050.viol*1000), mean(diff2099.viol*1000)) # change in mean daily violations per 1k inspections
meandiff.viol <- setZ(meandiff.viol, id) # set year as z variable
names(meandiff.viol) <- c('2050','2099') # set names
meandiff.viol <- trim(meandiff.viol, padding=0, values=NA) # trim outer NAs

# police projection----
theme <- rasterTheme(region=colorRampPalette(c("orange3","orange","gray80","cornflowerblue","steelblue","blue","blue4"))(100)) # define color palette

pol.proj <- levelplot(meandiff.pol, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               names=c('2050','2099'), # set names
               par.strip.text=list(cex=1.5), # set size of names
               layout=c(2,1), # horizontal layout
               colorkey=list(space="bottom", labels=list(cex=1.35)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
pol.proj <- update(pol.proj, par.settings = list(axis.line = list(col = 0)), sub=list("Projected Change in Mean Daily Police Stops\n(Percentage Points)", cex=1.5)) # remove plot borders and add subtitle below legend
pol.proj <- pol.proj + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# crash projection----
theme <- rasterTheme(region=colorRampPalette(c("gray90","lavenderblush","darksalmon","red1","red2","red3","red4"))(100)) # define color palette

cr.proj <- levelplot(meandiff.cr, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               names=c('2050','2099'), # set names
               par.strip.text=list(cex=1.5), # set size of names
               layout=c(2,1), # horizontal layout
               colorkey=list(space="bottom", labels=list(cex=1.35)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
cr.proj <- update(cr.proj, par.settings = list(axis.line = list(col = 0)), sub=list("Projected Change in Mean Daily Fatal Crash Risk\n(Percentage Points)", cex=1.5)) # remove plot borders and add subtitle below legend
cr.proj <- cr.proj + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# food safety inspection projection----
theme <- rasterTheme(region=colorRampPalette(c("darkorchid4","darkorchid","gray80","yellow","yellow2","gold","gold3"))(100)) # define color palette

insp.proj <- levelplot(meandiff.insp, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               names=c('2050','2099'), # set names
               par.strip.text=list(cex=1.5), # set size of names
               layout=c(2,1), # horizontal layout
               colorkey=list(space="bottom", labels=list(cex=1.35)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
insp.proj <- update(insp.proj, par.settings = list(axis.line = list(col = 0)), sub=list("Projected Change in Mean Daily Food Safety Inspections\n(per 100,000 Facilities)", cex=1.5)) # remove plot borders and add subtitle below legend
insp.proj <- insp.proj + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# food safety violation projection----
theme <- rasterTheme(region=colorRampPalette(c("gray95","darkolivegreen2","darkolivegreen3","darkolivegreen4","darkolivegreen"))(100)) # define color palette

viol.proj <- levelplot(meandiff.viol, 
               par.settings = theme, #set theme
               xlab=NULL, # kill x-axis
               ylab=NULL, # kill y-axis
               scales=list(draw=FALSE),
               names=c('2050','2099'), # set names
               par.strip.text=list(cex=1.5), # set size of names
               layout=c(2,1), # horizontal layout
               colorkey=list(space="bottom", labels=list(cex=1.35)), # make legend on bottom of plot
               panel=panel.levelplot.raster, # makes it so that grid lines are not plotted
               pretty=TRUE
    )
viol.proj <- update(viol.proj, par.settings = list(axis.line = list(col = 0)), sub=list("Projected Change in Mean Daily Food Safety Violations\n(per 1,000 Inspections)", cex=1.5)) # remove plot borders and add subtitle below legend
viol.proj <- viol.proj + layer(sp.lines(states, lwd=0.5, col="white")) # add us shapefiles

# plot----
cairo_pdf(filename = './figure4.pdf', width=14, height=7.25, fallback_resolution=1200)
grid.arrange(
    arrangeGrob(pol.proj, cr.proj, ncol=2),
    rectGrob(gp=gpar(col="white")),
    arrangeGrob(insp.proj, viol.proj, ncol=2), nrow=3, heights=c(0.48, 0.04, 0.48))
grid.text("a", x=unit(0.03, "npc"), y=unit(0.96, "npc"), gp=gpar(fontsize=25,font=2))
grid.text("b", x=unit(0.53, "npc"), y=unit(0.96, "npc"), gp=gpar(fontsize=25,font=2))
grid.text("c", x=unit(0.03, "npc"), y=unit(0.45, "npc"), gp=gpar(fontsize=25,font=2))
grid.text("d", x=unit(0.53, "npc"), y=unit(0.45, "npc"), gp=gpar(fontsize=25,font=2))
dev.off()
```

```{r warm, out.width="100%", fig.cap="Projected change in regulatory behaviors and outcomes in the US due to future warming. This figure presents the 25km x 25km grid cell projections of the effects of warming in the United States over this century on the regulatory outcomes examined in this study. Projections are calculated using downscaled climatic model data from NASA's NEX under the RCP8.5 high emissions scenario across the mean of the 21 CMIP5 models in the ensemble. We couple these climate model data with the estimates from our historical statistical models to project the mean effects of climate change on each outcome. Panel (a) depicts projected changes to police stops, panel (b) depicts projected changes to fatal crash risk, panel (c) depicts projected changes to food safety inspections, and panel (d) depicts projected changes to food safety violations. Annual averages smooth the seasonal variation in these outcomes (see SI: Climate impact projections for by-season projections)."}
knitr::include_graphics("figure4.pdf")
```

```{r}
warm <- fig_num(name = "warm", caption = "")
```
